function [parms_fitted, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
    FitModel20220926_comp_EU_stdU(parms, b_guess, kappa_guess, ...
        U_mean_target, U_std_target, NumTries, fsolve_opts)
% target moments:
% E[U], vol[u]
% parameters:
% b, kappa

    if b_guess<=0; error('check b_guess'); end
    if kappa_guess<=0; error('check kappa_guess'); end
    
    fixedpt_tol = 1e-6; % speed things up
    
    N_Nxz = repmat(parms.grid_N(:),[1 parms.Nx parms.Nz]);
    U_Nxz = 1 - N_Nxz;
    
    FLAG_SUCCESSFUL_FSOLVE = false;
    
    for iter = 1:NumTries
        fprintf('iter %g with b_guess=%g, kappa_guess=%g.\n', ...
            iter, b_guess, kappa_guess);
        [XFit, FitVal, EXITFLAG] = fsolve(@FitFun,[log(b_guess); log(kappa_guess)],fsolve_opts);
        if EXITFLAG>0
            fprintf('Found roots on iteration %g.\n', iter);
            FLAG_SUCCESSFUL_FSOLVE = true;
            break;
        else
            fprintf('iter %g, fitted vals (not solution):\n b=%g, kappa=%g.\n', ...
                iter, exp(XFit(1)), exp(XFit(2)));
            b_guess = 0.3 + 0.6*rand;
            kappa_guess = 2*rand;
        end
    end
    
    % do a last compute
    parms_fitted = parms;
    parms_fitted.b(:) = exp(XFit(1));    
    parms_fitted.kappa(:) = exp(XFit(2));

    fprintf('kappa=%g, b=%g, kappa=%g.\n', ...
        parms_fitted.b(1), parms_fitted.kappa(1));

    gesol = SolveFixedPoint_Model_20220926(parms_fitted,fixedpt_tol);
    parms_fitted.FitVal = FitVal;

    if ~FLAG_SUCCESSFUL_FSOLVE
        warning('Model fit: not all moments are satisfied.');
    end

    function Z = FitFun(X)
        
        Z = NaN*zeros(2,1);
        
        parms_temp = parms;
        parms_temp.b(:) = exp(X(1));
        parms_temp.kappa(:) = exp(X(2));

        gesol_temp = SolveFixedPoint_Model_20220926(parms_temp,fixedpt_tol);
        if gesol_temp.VFIerr>fixedpt_tol
            warning('FitModel20220926_ver1: VFI error = %g is greater than tolerance of %g.', ...
                gesol_temp.VFIerr, fixedpt_tol)
            bComputeMoments = false;
        else
            bComputeMoments = true;
        end
        
        %% Newer version
        
        if bComputeMoments

        [Phi_Nxz,ss_dist_Nxz] = MakeSparseTransMatrix_Nxz(parms_temp,gesol_temp);

        % NEW CODE: quarterly average
        [mean_U, std_U] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,U_Nxz);
        
        Z(1) = (mean_U - U_mean_target)/U_mean_target;
        Z(2) = (std_U - U_std_target)/U_std_target;
        
        end
        
    end

end